﻿/*
   Copyright 2010 University of Southampton

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
*/

/*
 * Computed Tomography voxel set handlers
 * 
 * Dr. Richard Boardman
 * Engineering Materials Group
 * School of Engineering Sciences
 * Copyright (c) 2010 University of Southampton
 *
 */

using System;
using System.IO;
using System.Windows.Media.Imaging;
using System.Windows.Media;
using System.Diagnostics;

namespace uk.ac.soton.ses
{
    /// <summary>
    /// General exception thrown by CtReader
    /// </summary>
    public class CtReaderException : Exception
    {
        /// <summary>
        /// A CtReader exception
        /// </summary>
        /// <param name="message">The message to throw</param>
        public CtReaderException(string message)
            : base(message)
        {
            // using parent constructor
        }
    }

    /// <summary>
    /// Provides handling methods for computed tomography data sets
    /// </summary>
    public class CtReader
    {
        /// <summary>
        /// Provides dimensionality information for voxel sets
        /// </summary>
        public class Dimensions
        {
            /// <summary>
            /// The x length of the voxel set
            /// </summary>
            public long xsize = 0;

            /// <summary>
            /// The y length of the voxel set
            /// </summary>
            public long ysize = 0;

            /// <summary>
            /// The z length of the voxel set
            /// </summary>
            public long zsize = 0;

            /// <summary>
            /// Creates a new dimension object with the provided sizes
            /// </summary>
            /// <param name="xsize">The X length of the voxel set</param>
            /// <param name="ysize">The Y length of the voxel set</param>
            /// <param name="zsize">The Z length of the voxel set</param>
            public Dimensions(int xsize, int ysize, int zsize)
            {
                this.xsize = xsize;
                this.ysize = ysize;
                this.zsize = zsize;
            }

            /// <summary>
            /// Gets the approximate size of the dataset in bytes
            /// </summary>
            /// <param name="depth">Dataset depth in bits</param>
            public long GetSize(int depth)
            {
                return this.xsize * this.ysize * this.zsize * (depth / 8);
            }

            /// <summary>
            /// Trivially compares the provided size with the dimensionality to compute the bit depth of the data
            /// </summary>
            /// <param name="size">The total size of the voxel set</param>
            /// <returns>Estimated bit depth</returns>
            public int GetBitDepth(long size)
            {
                for (int i = 0; i < 128; i++)
                {
                    if (size == this.GetSize(i))
                    {
                        return i;
                    }
                }
                return -1;
            }

            /// <summary>
            /// Gets the minimum uncompressed byte size of the voxel set based on the dimensions if the depth were 8-bit
            /// </summary>
            /// <returns>Voxel set size in bytes</returns>
            public long Get8BitSize()
            {
                return GetSize(8);
            }

            /// <summary>
            /// Gets the minimum uncompressed byte size of the voxel set based on the dimensions if the depth were 16-bit
            /// </summary>
            /// <returns>Voxel set size in bytes</returns>
            public long Get16BitSize()
            {
                return this.GetSize(16);
            }

            /// <summary>
            /// Gets the minimum uncompressed byte size of the voxel set based on the dimensions if the depth were 32-bit
            /// </summary>
            /// <returns>Voxel set size in bytes</returns>
            public long Get32BitSize()
            {
                return this.GetSize(32);
            }

            public long Get64BitSize()
            {
                return this.GetSize(64);
            }
        }

        /// <summary>
        /// Provides some colour look-up tables as BitmapPalettes
        /// </summary>
        public class LookupTableColours
        {
            // static methods

            /// <summary>
            /// Generates a red-to-blue palette
            /// </summary>
            /// <returns>A red-to-blue palette</returns>
            public static BitmapPalette RedBlue()
            {

                // range ... 255?
                Color[] c = new Color[256];

                for (int i = 0; i < 255; i++)
                {
                    c[i].R = (byte)(255 - i);
                    c[i].G = 20;
                    c[i].B = (byte)i;
                    c[i].A = 255;
                }                
                
                return new BitmapPalette(c);
            }

            
            /// <summary>
            /// Generates a green-to-purple palette
            /// </summary>
            /// <returns>A green-to-purple palette</returns>
            public static BitmapPalette GreenPurple()
            {

                // range ... 255?
                Color[] c = new Color[256];

                for (int i = 0; i < 255; i++)
                {
                    c[i].R = (byte)((255 - i) / 2);
                    c[i].G = (byte)i;
                    c[i].B = (byte)(255 - i);
                    c[i].A = 200;
                }

                return new BitmapPalette(c);
            }

            /// <summary>
            /// Generates a red-through-green-to-blue palette
            /// </summary>
            /// <returns>A red-through-green-to-blue palette</returns>
            public static BitmapPalette RedGreenBlue()
            {

                // range ... 255?
                Color[] c = new Color[256];

                for (int i = 0; i < 128; i++)
                {
                    c[i].R = (byte)(255 - (2 * i));
                    c[i].G = (byte)(2 * i);
                    c[i].B = 0;
                    c[i].A = 255;
                }

                for (int i = 128; i < 256; i++)
                {
                    c[i].R = 0;
                    c[i].G = (byte)(255 - (2 * i));
                    c[i].B = (byte)(2 * i);
                    c[i].A = 255;
                }

                return new BitmapPalette(c);
            }

            /// <summary>
            /// Generates a random colour palette
            /// </summary>
            /// <returns>A random colour palette</returns>
            public static BitmapPalette Random()
            {

                Random r = new Random();
                // range ... 255?
                Color[] c = new Color[256];

                for (int i = 0; i < 256; i++)
                {
                    c[i].R = (byte)(r.Next(256));
                    c[i].G = (byte)(r.Next(256));
                    c[i].B = (byte)(r.Next(256));
                    c[i].A = 255; // (byte)(r.Next(256));
                }

                return new BitmapPalette(c);
            }

        }

        /// <summary>
        /// Planes for CT data slices: XY, XZ or YZ
        /// </summary>
        public enum Plane
        {
            /// <summary>
            /// The X-Y plane, with slices in Z
            /// </summary>
            XY,
            /// <summary>
            /// The X-Z plane, with slices in Y
            /// </summary>
            XZ,
            /// <summary>
            /// The Y-Z plane, with slices in X
            /// </summary>
            YZ
        }

        /// <summary>
        /// The last call to the Message method
        /// </summary>
        private DateTime lastMessageDateTime = DateTime.Now;

        /// <summary>
        /// The raw voxel data source
        /// </summary>
        public string sourceFile = null;

        /// <summary>
        /// The size of the voxel data source
        /// </summary>
        private long fileSize = -1;

        /// <summary>
        /// The dimensions of the voxel data source
        /// </summary>
        private Dimensions dimensions = null;

        /// <summary>
        /// Creates a new CtReader from the raw voxel data in sourceFile
        /// </summary>
        /// <param name="sourceFile">Raw voxel data source file</param>
        public CtReader(string sourceFile, Dimensions dimensions)
        {
            this.sourceFile = sourceFile;
            this.dimensions = dimensions;
            this.fileSize = CalculateFileSize();
        }

        /// <summary>
        /// The output DPI of any generated BitmapSources or images
        /// </summary>
        public static int OUTPUT_IMAGE_DPI = 96;

        /// <summary>
        /// Gets the bit depth of the voxel values
        /// </summary>
        /// <returns>Bit depth</returns>
        public int GetDepth()
        {
            if (this.dimensions != null)
            {
                return this.dimensions.GetBitDepth(this.FileSize);
            }
            throw new CtReaderException("The dimensions have not been set for this file");
        }

        /// <summary>
        /// Returns the X size of the associated dimension
        /// </summary>
        /// <returns>X size</returns>
        public long GetXSize()
        {
            return this.dimensions.xsize;
        }

        /// <summary>
        /// Returns the Y size of the associated dimension
        /// </summary>
        /// <returns>Y size</returns>
        public long GetYSize()
        {
            return this.dimensions.ysize;
        }

        /// <summary>
        /// Returns the Z size of the associated dimension
        /// </summary>
        /// <returns>Z size</returns>
        public long GetZSize()
        {
            return this.dimensions.zsize;
        }

        /// <summary>
        /// Returns a <code>CtReaderStream</code> representing the current instance
        /// </summary>
        /// <returns><code>CtReaderStream</code></returns>
        public CtReaderStream GetReaderStream()
        {            
           return new CtReaderStream(this);
        }

        /// <summary>
        /// Returns a <code>CtReaderStream</code> representing the current instance, clipped to the cube defined by the provided parameters
        /// </summary>
        /// <param name="low_x">The low X value for clipping</param>
        /// <param name="low_y">The low Y value for clipping</param>
        /// <param name="low_z">The low Z value for clipping</param>
        /// <param name="high_x">The high X value for clipping</param>
        /// <param name="high_y">The high Y value for clipping</param>
        /// <param name="high_z">The high Z value for clipping</param>
        /// <returns><code>CtReaderStream</code></returns>
        public CtReaderStream GetReaderStream(long low_x, long low_y, long low_z, long high_x, long high_y, long high_z)
        {            
            return new CtReaderStream(this, low_x, low_y, low_z, high_x, high_y, high_z);
        }

        /// <summary>
        /// Sets the size of the voxel data source from the file data
        /// </summary>
        private long CalculateFileSize()
        {
            using (BinaryReader br = new BinaryReader(File.OpenRead(this.sourceFile)))
            {
                return br.BaseStream.Length;
            }
        }

        /// <summary>
        /// Returns the size of the voxel data source
        /// </summary>
        /// <returns>File size in bytes</returns>
        public long FileSize
        {
            get
            {
                if (fileSize != -1)
                {
                    return this.fileSize;
                }
                else
                    throw new CtReaderException("File size not set");
            }
        }

        /// <summary>
        /// Puts a message to the standard output, with a timestamp and a time difference between this and the previous message
        /// </summary>
        /// <param name="message">The message to output</param>
        [Obsolete]
        public void Message(string message)
        {
            Trace.WriteLine(String.Format("[{0}] ({1}): {2}", DateTime.Now, DateTime.Now - this.lastMessageDateTime, message));
            this.lastMessageDateTime = DateTime.Now;
        }

        /// <summary>
        /// Gets the minimum value in the voxel data set (8-bit). Warning: this is a slow (unoptimised) method
        /// </summary>
        /// <returns>Minimum value</returns>
        public byte GetMinValue8()
        {
            byte b = Byte.MaxValue;
            byte t = b;
            BinaryReader br = null;
            long counter = 0;

            using (br = new BinaryReader(File.OpenRead(this.sourceFile)))
            {
                while (counter < this.fileSize)
                {
                    counter++;
                    t = br.ReadByte();
                    if (t < b)
                    {
                        b = t;
                    }
                }
            }            
            return b;
        }

        /// <summary>
        /// Gets the minimum value in the voxel data set (16-bit). Warning: this is a slow (unoptimised) method
        /// </summary>
        /// <returns>Minimum value</returns>
        public ushort GetMinValue16()
        {
            ushort b = UInt16.MaxValue;
            ushort t = b;
            BinaryReader br = null;
            long counter = 0;

            using (br = new BinaryReader(File.OpenRead(this.sourceFile)))
            {
                while (counter < this.fileSize)
                {
                    counter+=2;
                    t = br.ReadUInt16();
                    if (t < b)
                    {
                        b = t;
                    }
                }
            }            
            return b;
        }

        /// <summary>
        /// Gets the minimum value in the voxel data set (32-bit). Warning: this is a slow (unoptimised) method
        /// </summary>
        /// <returns>Minimum value</returns>
        public uint GetMinValue32()
        {
            uint b = UInt32.MaxValue;
            uint t = b;
            BinaryReader br = null;
            long counter = 0;

            using (br = new BinaryReader(File.OpenRead(this.sourceFile)))
            {
                while (counter < this.fileSize)
                {
                    counter+=4;
                    t = br.ReadUInt32();
                    if (t < b)
                    {
                        b = t;
                    }
                }
            }            
            return b;
        }
        /// <summary>
        /// Gets the maximum value in the voxel data set (8-bit). Warning: this is a slow (unoptimised) method
        /// </summary>
        /// <returns>Maximum value</returns>
        public byte GetMaxValue8()
        {
            byte b = Byte.MinValue;
            byte t = b;
            BinaryReader br = null;
            long counter = 0;
            using (br = new BinaryReader(File.OpenRead(this.sourceFile)))
            {
                while (counter < this.fileSize)
                {
                    counter++;
                    t = br.ReadByte();
                    if (t > b)
                    {
                        b = t;
                    }
                }
            }
            return b;
        }

        /// <summary>
        /// Gets the maximum value in the voxel data set (16-bit). Warning: this is a slow (unoptimised) method
        /// </summary>
        /// <returns>Maximum value</returns>
        public ushort GetMaxValue16()
        {
            ushort b = UInt16.MinValue;
            ushort t = b;
            BinaryReader br = null;
            long counter = 0;
            using (br = new BinaryReader(File.OpenRead(this.sourceFile)))
            {
                while (counter < this.fileSize)
                {
                    counter+=2;
                    t = br.ReadUInt16();
                    if (t > b)
                    {
                        b = t;
                    }
                }
            }            
            return b;
        }

        /// <summary>
        /// Gets the maximum value in the voxel data set (32-bit). Warning: this is a slow (unoptimised) method
        /// </summary>
        /// <returns>Maximum value</returns>
        public uint GetMaxValue32()
        {
            uint b = UInt32.MinValue;
            uint t = b;
            BinaryReader br = null;
            long counter = 0;
            using (br = new BinaryReader(File.OpenRead(this.sourceFile)))
            {
                while (counter < this.fileSize)
                {
                    counter+=4;
                    t = br.ReadUInt32();
                    if (t > b)
                    {
                        b = t;
                    }
                }
            }            
            return b;
        }

        public double GetMinValueDouble()
        {
            double b = Double.MaxValue;
            double t = b;
            BinaryReader br = null;
            long counter = 0;
            using (br = new BinaryReader(File.OpenRead(this.sourceFile)))
            {
                while (counter < this.fileSize)
                {
                    counter += 8;
                    t = br.ReadDouble();
                    if (t < b)
                    {
                        b = t;
                    }
                }
            }
            return b;
        }


        public double GetMaxValueDouble()
        {
            double b = Double.MinValue;
            double t = b;
            BinaryReader br = null;
            long counter = 0;
            using (br = new BinaryReader(File.OpenRead(this.sourceFile)))
            {
                while (counter < this.fileSize)
                {
                    counter += 8;
                    t = br.ReadDouble();
                    if (t > b)
                    {
                        b = t;
                    }
                }
            }
            return b;
        }

        /// <summary>
        /// Automatically scales the plane between 0 and 255 in the case of a byte plane. 
        /// 
        /// If we wanted to scale the data between 0 and 10 for the sake of example, and the plane contains
        /// the values:
        /// 
        /// <code>3 7 8</code>
        /// 
        /// then the lowest value (3) is removed from every value, and scaled upwards so the largest value is 10 (i.e.,
        /// multiplied by 10/(max-min) -> 10/(8-3) -> 10/5 -> 2. The resulting dataset would be:
        /// 
        /// <code>0 8 10</code>
        /// 
        /// </summary>
        /// <param name="plane">The byte array to scale</param>
        /// <returns>Scaled byte array</returns>
        public byte[] AutoContrastPlane(byte[] plane)
        {
            //byte[] contrastplane = new byte[plane.Length];

            byte[] interimplane = new byte[plane.Length];

            // max and min
            byte min = Byte.MaxValue;
            byte max = Byte.MinValue;

            for (int i = 0; i < plane.Length; i++)
            {
                if (plane[i] < min)
                {
                    min = plane[i];
                }
                if (plane[i] > max)
                {
                    max = plane[i];
                }
            }

            double scalefactor = ((double)Byte.MaxValue) / (max - min);
            double temp = 0;

            for (int i = 0; i < plane.Length; i++)
            {
                temp = ((((double)plane[i]) - min) * scalefactor);

                if (temp < 0)
                {
                    temp = 0;
                }
                if (temp > Byte.MaxValue)
                {
                    temp = Byte.MaxValue;
                }

                interimplane[i] = (byte)temp;
            }

            return interimplane;
        }

        /// <summary>
        /// Automatically scales the plane between 0 and 65535 in the case of a short plane. 
        /// </summary>
        /// <param name="plane">The short array to scale</param>
        /// <returns>Scaled short array</returns>
        public ushort[] AutoContrastPlane(ushort[] plane)
        {
            ushort[] interimplane = new ushort[plane.Length];

            // max and min
            ushort min = UInt16.MaxValue;
            ushort max = UInt16.MinValue;

            for (int i = 0; i < plane.Length; i++)
            {
                if (plane[i] < min)
                {
                    min = plane[i];
                }
                if (plane[i] > max)
                {
                    max = plane[i];
                }
            }

            double scalefactor = ((double)Int16.MaxValue) / (max - min);
            double temp = 0;

            for (int i = 0; i < plane.Length; i++)
            {
                temp = ((((double)plane[i]) - min) * scalefactor);

                if (temp < 0)
                {
                    temp = 0;
                }
                if (temp > UInt16.MaxValue)
                {
                    temp = UInt16.MaxValue;
                }

                interimplane[i] = (ushort)temp;
            }

            return interimplane;
        }

        /// <summary>
        /// Automatically scales the plane between 0 and 4294967295 in the case of an integer plane. 
        /// </summary>
        /// <param name="plane">The integer array to scale</param>
        /// <returns>Scaled integer array</returns>
        public uint[] AutoContrastPlane(uint[] plane)
        {
            uint[] interimplane = new uint[plane.Length];

            // max and min
            uint min = UInt32.MaxValue;
            uint max = UInt32.MaxValue;

            for (int i = 0; i < plane.Length; i++)
            {
                if (plane[i] < min)
                {
                    min = plane[i];
                }
                if (plane[i] > max)
                {
                    max = plane[i];
                }
            }

            double scalefactor = ((double)Int16.MaxValue) / (max - min);
            double temp = 0;

            for (int i = 0; i < plane.Length; i++)
            {
                temp = ((((double)plane[i]) - min) * scalefactor);

                if (temp < 0)
                {
                    temp = 0;
                }
                if (temp > Int16.MaxValue)
                {
                    temp = Int16.MaxValue;
                }

                interimplane[i] = (uint)temp;
            }

            return interimplane;
        }

        /// <summary>
        /// Gets the number of slices of data available on the specified plane
        /// </summary>
        /// <param name="planetype">The plane to consider</param>
        /// <returns></returns>
        public int GetSliceCount(CtReader.Plane planetype)
        {
            switch (planetype)
            {
                case Plane.XY:
                    return (int)this.dimensions.zsize;
                case Plane.XZ:
                    return (int)this.dimensions.ysize;
                case Plane.YZ:
                    return (int)this.dimensions.xsize;
                default:
                    throw new CtReaderException("Invalid plane supplied");
            }
        }

        /// <summary>
        /// Extracts a slice of interest across the specified plane from the voxel set. Slices are zero-indexed.
        /// </summary>
        /// <param name="slice">The index of the slice to extract</param>
        /// <param name="planetype">The plane type (XY, XZ or YZ) to extract</param>
        /// <param name="autoContrast">Whether the <code>AutoContrastPlane</code> method should be called</param>
        /// <returns>Slice</returns>
        public byte[] ExtractPlane(int slice, CtReader.Plane planetype, bool autoContrast)
        {
            return this.ExtractBytePlane(slice, planetype, autoContrast);
        }

        public uint[] ExtractPlane32(int slice, CtReader.Plane planetype, bool autoContrast)
        {
            if (autoContrast)
            {
                return this.AutoContrastPlane(this.ExtractArbitraryPlane(slice, planetype, 32));
            }
            return this.ExtractArbitraryPlane(slice, planetype, 32);
        }

        private ushort[] convertToUShort(uint[] data)
        {
            ushort[] converted = new ushort[data.Length];
            for (long l = 0; l < data.Length; l++)
            {
                if (data[l] > UInt16.MaxValue)
                {
                    throw new CtReaderException("Data range out of bounds for unsigned short");
                }
                converted[l] = (ushort)data[l];
            }
            return converted;
        }

        public ushort[] ExtractPlane16(int slice, CtReader.Plane planetype, bool autoContrast)
        {
            uint[] thisplane;

            if (autoContrast)
            {
                thisplane = this.AutoContrastPlane(this.ExtractArbitraryPlane(slice, planetype, 16));
            }
            else
            {
                thisplane = this.ExtractArbitraryPlane(slice, planetype, 16);
            }

            return this.convertToUShort(thisplane);
        }

        /// <summary>
        /// *** NOT YET TESTED ***
        /// Extended method to handle 16- and 32-bit voxel sets. These are *unsigned* at the moment but this may
        /// be adjusted to signed (or signing handling added) depending on the CT output.
        /// 
        /// Extracts a slice of interest across the specified plane from the voxel set. Slices are zero-indexed.
        /// </summary>
        /// <param name="slice">The index of the slice to extrac</param>
        /// <param name="planetype">The plane type (XY, XZ or YZ) to extract</param>
        /// <param name="autoContrast"></param>
        /// <param name="bitsize">The bit depth of the voxel set (16 or 32)</param>
        /// <returns></returns>
        public uint[] ExtractArbitraryPlane(int slice, CtReader.Plane planetype, int bitsize)
        {
            uint[] planedata;

            int seekstep = 1;

            switch (bitsize)
            {
                case 8:
                    throw new CtReaderException("Please use the ExtractBytePlane method");                    
                case 16:
                    seekstep = 2;
                    break;
                case 32:
                    seekstep = 4;
                    break;
                default:
                    throw new CtReaderException(String.Format("Cannot extract plane size of {0} bits; 8-, 16- and 32-bit are supported"));
            }

            using (BinaryReader br = new BinaryReader(File.OpenRead(this.sourceFile)))
            {
                switch (planetype)
                {
                    case CtReader.Plane.XY:
                        planedata = new uint[this.dimensions.xsize * this.dimensions.ysize];
                        br.BaseStream.Seek(this.dimensions.xsize * this.dimensions.ysize * slice * seekstep, SeekOrigin.Begin);
                        for (long i = 0; i < planedata.Length; i++)
                        {
                            if (bitsize == 16)
                            {
                                planedata[i] = (uint)br.ReadUInt16();
                            }
                            else if (bitsize == 32)
                            {
                                planedata[i] = br.ReadUInt32();                                
                            }
                        }                        
                        break;
                    case CtReader.Plane.XZ:
                        planedata = new uint[this.dimensions.xsize * this.dimensions.zsize];
                        for (long i = 0; i < this.dimensions.zsize; i++)
                        {
                            br.BaseStream.Seek((slice * this.dimensions.xsize * seekstep) + (this.dimensions.xsize * this.dimensions.ysize * i * seekstep), SeekOrigin.Begin);
                            for (long j = 0; j < this.dimensions.xsize; j++)
                            {
                                if (bitsize == 16)
                                {
                                    planedata[j + (i * this.dimensions.xsize)] = br.ReadUInt16();
                                }
                                else if (bitsize == 32)
                                {
                                    planedata[j + (i * this.dimensions.xsize)] = br.ReadUInt32();
                                }
                            }
                        }
                        break;                        
                    case CtReader.Plane.YZ:
                        planedata = new uint[this.dimensions.ysize * this.dimensions.zsize];
                                                
                        for (long i = 0; i < this.dimensions.zsize; i++)
                        {
                            for (long j = 0; j < this.dimensions.ysize; j++)
                            {
                                br.BaseStream.Seek((seekstep * slice) +
                                                  (seekstep * this.dimensions.xsize * j) +
                                                  (seekstep * this.dimensions.xsize * this.dimensions.ysize * i),
                                                   SeekOrigin.Begin);
                                if (bitsize == 16)
                                {
                                    planedata[j + (i * this.dimensions.ysize)] = br.ReadUInt16();
                                }
                                else if (bitsize == 32)
                                {
                                    planedata[j + (i * this.dimensions.ysize)] = br.ReadUInt32();
                                }
                            }
                        }
                        break;
                    default:
                        throw new CtReaderException("Unknown plane type");
                }
            }

            return planedata;
        }

        public double[] ExtractFloatPlane(int slice, CtReader.Plane planetype, int bitsize)
        {
            double[] planedata;

            int seekstep = 1;

            switch (bitsize)
            {
                case 8:
                    throw new CtReaderException("Please use the ExtractBytePlane method");
                case 16:                    
                    throw new CtReaderException("16-bit (2-byte) floating point numbers are not supported");
                case 32:
                    seekstep = 4;
                    break;
                case 64:
                    seekstep = 8;
                    break;
                default:
                    throw new CtReaderException(String.Format("Cannot extract plane size of {0} bits; 8-, 16- and 32-bit are supported"));
            }

            using (BinaryReader br = new BinaryReader(File.OpenRead(this.sourceFile)))
            {
                switch (planetype)
                {
                    case CtReader.Plane.XY:
                        planedata = new double[this.dimensions.xsize * this.dimensions.ysize];
                        br.BaseStream.Seek(this.dimensions.xsize * this.dimensions.ysize * slice * seekstep, SeekOrigin.Begin);
                        for (long i = 0; i < planedata.Length; i++)
                        {
                            if (bitsize == 32)
                            {
                                try
                                {
                                    planedata[i] = br.ReadSingle();
                                }
                                catch
                                {
                                    planedata[i] = 0.0;
                                    break;
                                }
                            }
                            else if (bitsize == 64)
                            {
                                planedata[i] = br.ReadDouble();
                            }
                        }
                        break;
                    case CtReader.Plane.XZ:
                        planedata = new double[this.dimensions.xsize * this.dimensions.zsize];
                        for (long i = 0; i < this.dimensions.zsize; i++)
                        {
                            br.BaseStream.Seek((slice * this.dimensions.xsize * seekstep) + (this.dimensions.xsize * this.dimensions.ysize * i * seekstep), SeekOrigin.Begin);
                            for (long j = 0; j < this.dimensions.xsize; j++)
                            {
                                if (bitsize == 32)
                                {
                                    planedata[j + (i * this.dimensions.xsize)] = (double)br.ReadSingle();
                                }
                                else if (bitsize == 64)
                                {
                                    planedata[j + (i * this.dimensions.xsize)] = br.ReadDouble();
                                }
                            }
                        }
                        break;
                    case CtReader.Plane.YZ:
                        planedata = new double[this.dimensions.ysize * this.dimensions.zsize];

                        for (long i = 0; i < this.dimensions.zsize; i++)
                        {
                            for (long j = 0; j < this.dimensions.ysize; j++)
                            {
                                br.BaseStream.Seek((seekstep * slice) +
                                                  (seekstep * this.dimensions.xsize * j) +
                                                  (seekstep * this.dimensions.xsize * this.dimensions.ysize * i),
                                                   SeekOrigin.Begin);
                                if (bitsize == 32)
                                {
                                    planedata[j + (i * this.dimensions.ysize)] = (double)br.ReadSingle();
                                }
                                else if (bitsize == 64)
                                {
                                    planedata[j + (i * this.dimensions.ysize)] = br.ReadDouble();
                                }
                            }
                        }
                        break;
                    default:
                        throw new CtReaderException("Unknown plane type");
                }
            }

            return planedata;
        }

        /// <summary>
        /// Extracts a slice of interest across the specified plane from the voxel set. Slices are zero-indexed.
        /// </summary>
        /// <param name="slice">The index of the slice to extract</param>
        /// <param name="planetype">The plane type (XY, XZ or YZ) to extract</param>
        /// <returns>Slice</returns>
        public byte[] ExtractBytePlane(int slice, CtReader.Plane planetype, bool autoContrast)
        {         
            byte[] planebytes;

            // open the stream
            BinaryReader br = null;
            using (br = new BinaryReader(File.OpenRead(this.sourceFile)))
            {
                switch (planetype)
                {
                    case CtReader.Plane.XY:
                        planebytes = new byte[this.dimensions.xsize * this.dimensions.ysize];
                        br.BaseStream.Seek(this.dimensions.xsize * this.dimensions.ysize * slice, SeekOrigin.Begin);
                        planebytes = br.ReadBytes((int)this.dimensions.xsize * (int)this.dimensions.ysize);                        
                        break;
                    case CtReader.Plane.XZ:
                        planebytes = new byte[this.dimensions.xsize * this.dimensions.zsize];
                        for (long i = 0; i < this.dimensions.zsize; i++)
                        {
                            br.BaseStream.Seek((slice * this.dimensions.xsize) + (this.dimensions.xsize * this.dimensions.ysize * i), SeekOrigin.Begin);
                            for (long j = 0; j < this.dimensions.xsize; j++)
                            {
                                planebytes[j + (i * this.dimensions.xsize)] = br.ReadByte();
                            }
                        }
                        break;
                    case CtReader.Plane.YZ:
                        planebytes = new byte[this.dimensions.ysize * this.dimensions.zsize];
                        for (long i = 0; i < this.dimensions.zsize; i++)
                        {
                            for (long j = 0; j < this.dimensions.ysize; j++)
                            {
                                br.BaseStream.Seek(slice +
                                                  (this.dimensions.xsize * j) +
                                                  (this.dimensions.xsize * this.dimensions.ysize * i),
                                                   SeekOrigin.Begin);
                                planebytes[j + (i * this.dimensions.ysize)] = br.ReadByte();
                            }
                        }
                        break;
                    default:
                        throw new CtReaderException("Unknown plane type");
                }
            }

            return (autoContrast ? this.AutoContrastPlane(planebytes) : planebytes);
        }

        /// <summary>
        /// Gets the correct X size for an output image on the proposed plane
        /// </summary>
        /// <param name="planetype">The plane type</param>
        /// <returns>Target bitmap X size</returns>
        private int GetPlaneXSize(CtReader.Plane planetype)
        {
            if (planetype == CtReader.Plane.XY || planetype == CtReader.Plane.XZ)
            {
                return (int)this.dimensions.xsize;
            }
            if (planetype == CtReader.Plane.YZ)
            {
                return (int)this.dimensions.ysize;
            }
            throw new CtReaderException("Unknown plane type");
        }

        /// <summary>
        /// Gets the correct Y size for an output image on the proposed plane
        /// </summary>
        /// <param name="planetype">The plane type</param>
        /// <returns>Target bitmap Y size</returns>
        private int GetPlaneYSize(CtReader.Plane planetype)
        {
            if (planetype == CtReader.Plane.XZ || planetype == CtReader.Plane.YZ)
            {
                return (int)this.dimensions.zsize;
            }
            if (planetype == CtReader.Plane.XY)
            {
                return (int)this.dimensions.ysize;
            }
            throw new CtReaderException("Unknown plane type");
        }

        /// <summary>
        /// Writes a TIFF file to filename from the supplied BitmapSource
        /// </summary>
        /// <param name="bs">The BitmapSource to use</param>
        /// <param name="filename">The filename to write</param>
        public void WriteTiff(BitmapSource bs, string filename)
        {            
            TiffBitmapEncoder tbe = new TiffBitmapEncoder();
            tbe.Frames.Add(BitmapFrame.Create(bs));
            tbe.Save(File.Open(filename, FileMode.Create));
        }

        /// <summary>
        /// Writes a JPEG file to filename from the supplied BitmapSource
        /// </summary>
        /// <param name="bs">The BitmapSource to use</param>
        /// <param name="filename">The filename to write</param>
        public void WriteJpeg(BitmapSource bs, string filename)
        {            
            JpegBitmapEncoder jbe = new JpegBitmapEncoder();
            jbe.Frames.Add(BitmapFrame.Create(bs));
            jbe.Save(File.Open(filename, FileMode.Create));
        }

        /// <summary>
        /// Writes a PNG file to filename from the supplied BitmapSource
        /// </summary>
        /// <param name="bs">The BitmapSource to use</param>
        /// <param name="filename">The filename to write</param>
        public void WritePng(BitmapSource bs, string filename)
        {            
            PngBitmapEncoder pbe = new PngBitmapEncoder();
            pbe.Frames.Add(BitmapFrame.Create(bs));
            pbe.Save(File.Open(filename, FileMode.Create));

        }

        /// <summary>
        /// Gets the plane slice provided as a BitmapSource appropriate for the supplied plane type. 
        /// Defaults to a greyscale palette.
        /// </summary>
        /// <param name="plane">The plane slice data</param>
        /// <param name="planetype">The plane type</param>
        /// <returns>A bitmap source</returns>
        public BitmapSource GetPlaneAsBitmap(byte[] plane, CtReader.Plane planetype)
        {
            //BitmapPalette bp = LookupTableColours.GreenPurple();
            BitmapPalette bp = BitmapPalettes.Gray256; 
            return this.GetPlaneAsBitmap(plane, planetype, bp);
        }

        /// <summary>
        /// Gets the plane slice provided as a BitmapSource appropriate for the supplied plane type.
        /// *** WARNING *** the palette cannot be defined in 32bpp (128-bit colour) so a <code>null</code>
        /// is used here. Behaviour may be undefined
        /// </summary>
        /// <param name="plane">The plane slice data</param>
        /// <param name="planetype">The plane type</param>
        /// <returns>A bitmap source</returns>
        public BitmapSource GetPlaneAsBitmap(ushort[] plane, CtReader.Plane planetype)
        {
            return this.GetPlaneAsBitmap(plane, planetype, null);
        }

        /// <summary>
        /// Gets the plane slice provided as a BitmapSource appropriate for the supplied plane type.
        /// *** WARNING *** the palette cannot be defined in 32bpp (128-bit colour) so a <code>null</code>
        /// is used here. Behaviour may be undefined
        /// </summary>
        /// <param name="plane">The plane slice data</param>
        /// <param name="planetype">The plane type</param>
        /// <returns>A bitmap source</returns>
        public BitmapSource GetPlaneAsBitmap(uint[] plane, CtReader.Plane planetype)
        {
            return this.GetPlaneAsBitmap(plane, planetype, null);
        }

        /// <summary>
        /// Gets the plane slice provided as a BitmapSource appropriate for the supplied plane type.
        /// </summary>
        /// <param name="plane">The plane slice data</param>
        /// <param name="planetype">The plane type</param>
        /// <param name="bp">The palette to use</param>
        public BitmapSource GetPlaneAsBitmap(byte[] plane, CtReader.Plane planetype, BitmapPalette bp)
        {            
            int xsize = this.GetPlaneXSize(planetype);
            int ysize = this.GetPlaneYSize(planetype);

            BitmapSource bs = BitmapSource.Create(xsize,
                                                  ysize,
                                                  OUTPUT_IMAGE_DPI,
                                                  OUTPUT_IMAGE_DPI,
                                                  PixelFormats.Indexed8,
                                                  bp,
                                                  plane,
                                                  xsize);

            return bs;

        }

       

        /// <summary>
        /// Gets the plane slice provided as a BitmapSource appropriate for the supplied plane type.
        /// </summary>
        /// <param name="plane">The plane slice data</param>
        /// <param name="planetype">The plane type</param>
        /// <param name="bp">The palette to use</param>
        public BitmapSource GetPlaneAsBitmap(ushort[] plane, CtReader.Plane planetype, BitmapPalette bp)
        {                        
            int xsize = this.GetPlaneXSize(planetype);
            int ysize = this.GetPlaneYSize(planetype);
            return BitmapSource.Create(xsize,
                                       ysize,
                                       OUTPUT_IMAGE_DPI,
                                       OUTPUT_IMAGE_DPI,
                                       PixelFormats.Gray16,
                                       bp,
                                       plane,
                                       xsize*2);            
        }

        /// <summary>
        /// Gets the plane slice provided as a BitmapSource appropriate for the supplied plane type.
        /// </summary>
        /// <param name="plane">The plane slice data</param>
        /// <param name="planetype">The plane type</param>
        /// <param name="bp">The palette to use</param>
        public BitmapSource GetPlaneAsBitmap(uint[] plane, CtReader.Plane planetype, BitmapPalette bp)
        {            
            int xsize = this.GetPlaneXSize(planetype);
            int ysize = this.GetPlaneYSize(planetype); 
           
            // convert plane to float?
            double[] floatplane = new double[plane.LongLength];
            for (long l = 0; l < plane.LongLength; l++)
            {
               floatplane[l] = (double)plane[l] / (double)uint.MaxValue;
            }

            BitmapSource bs = BitmapSource.Create(xsize,
                                                  ysize,
                                                  OUTPUT_IMAGE_DPI,
                                                  OUTPUT_IMAGE_DPI,                                                  
                                                  PixelFormats.Gray32Float,
                                                  bp,
                                                  floatplane,
                                                  xsize*4);

            return bs;
        }

        /// <summary>
        /// Casts the incoming double array as a float array on an element-wise basis. Note that this may
        /// result in a loss of precision.
        /// </summary>
        /// <param name="arrayin">The double array to cast as float</param>
        /// <returns>Float array</returns>
        private static float[] CastDoubleArrayAsFloatArray(double[] arrayin)
        {
            float[] arrayout = new float[arrayin.LongLength];
            for (long l = 0; l < arrayin.LongLength; l++)
            {
                arrayout[l] = (float)arrayin[l];
            }            
            return arrayout;
        }

        /// <summary>
        /// Gets the plane slice provided as a BitmapSource appropriate for the supplied plane type.
        /// </summary>
        /// <param name="plane">Double array containing the slice information</param>
        /// <param name="planetype">The plane type</param>
        /// <param name="bp">The BitmapPalette to use (<code>null</code> for default)</param>
        /// <param name="bitsize">The bit size of the incoming data (64 for double, 32 for float); if 32 (float) is requested the supplied double array will be downcast</param>
        /// <param name="useAutoContrasting"><code>true</code> to enable autocontrasting, else <code>false</code> - <code>true</code> is recommended for imaging unless the bounds
        /// for scaling are known</param>        
        /// <returns>BitmapSource representing the plane data</returns>
        public BitmapSource GetPlaneAsBitmap(double[] plane, CtReader.Plane planetype, BitmapPalette bp, int bitsize, bool useAutoContrasting)
        {
            switch (bitsize)
            {
                case 32:
                    return GetPlaneAsBitmap(CastDoubleArrayAsFloatArray(plane), planetype, bp, useAutoContrasting);
                case 64:
                    return GetPlaneAsBitmap(plane, planetype, bp, useAutoContrasting);
                default:
                    throw new CtReaderException(String.Format("Cannot use bit depth of {0} with floating-point planes; use 32- or 64-bit instead", bitsize));
            }
        }

        /// <summary>
        /// Gets a BitmapSource representing the supplied float[] array
        /// </summary>
        /// <param name="plane">The plane data</param>
        /// <param name="planetype">The plane type</param>
        /// <param name="bp">The BitmapPalette to use; <code>null</code> to use the default</param>
        /// <param name="useAutoContrasting"><code>true</code> to enable autocontrasting, else <code>false</code> - <code>true</code> is recommended for imaging unless the bounds
        /// for scaling are known</param>
        /// <returns>BitmapSource representing the plane data</returns>
        public BitmapSource GetPlaneAsBitmap(float[] plane, CtReader.Plane planetype, BitmapPalette bp, bool useAutoContrasting)
        {            
            int xsize = this.GetPlaneXSize(planetype);
            int ysize = this.GetPlaneYSize(planetype);
            int stridefactor = 4;
            float fmax = float.MinValue;
            float fmin = float.MaxValue;

            if (useAutoContrasting)
            {
                for (long l = 0; l < plane.LongLength; l++)
                {
                    if (plane[l] > fmax)
                    {
                        fmax = plane[l];
                    }
                    if (plane[l] < fmin)
                    {
                        fmin = plane[l];
                    }
                }

                float foff = fmax + fmin;

                for (long l = 0; l < plane.LongLength; l++)
                {
                    plane[l] = (plane[l] + fmin) / foff;                    
                }
            }
            return BitmapSource.Create(xsize,
                                       ysize,
                                       OUTPUT_IMAGE_DPI,
                                       OUTPUT_IMAGE_DPI,
                                       PixelFormats.Gray32Float,
                                       bp,
                                       plane,
                                       xsize * stridefactor);
        }

        /// <summary>
        /// Gets a BitmapSource representing the supplied double[] array
        /// </summary>
        /// <param name="plane">The plane data</param>
        /// <param name="planetype">The plane type</param>
        /// <param name="bp">The BitmapPalette to use; <code>null</code> to use the default</param>
        /// <param name="useAutoContrasting"><code>true</code> to enable autocontrasting, else <code>false</code> - <code>true</code> is recommended for imaging unless the bounds
        /// for scaling are known</param>
        /// <returns>BitmapSource representing the plane data</returns>
        public BitmapSource GetPlaneAsBitmap(double[] plane, CtReader.Plane planetype, BitmapPalette bp, bool useAutoContrasting)
        {
            int xsize = this.GetPlaneXSize(planetype);
            int ysize = this.GetPlaneYSize(planetype);            
            int stridefactor = 8;
            
            if (useAutoContrasting)
            {
                double vmax = Double.MinValue;
                double vmin = Double.MaxValue;

                for (long l = 0; l < plane.LongLength; l++)
                {
                    if (plane[l] > vmax)
                    {
                        vmax = plane[l];
                    }
                    if (plane[l] < vmin)
                    {
                        vmin = plane[l];
                    }
                }

                double voff = vmin + vmax;

                for (long l = 0; l < plane.LongLength; l++)
                {                  
                    plane[l] = (plane[l] + vmin) / voff;
                }                
            }

            return BitmapSource.Create(xsize,
                                      ysize,
                                      OUTPUT_IMAGE_DPI,
                                      OUTPUT_IMAGE_DPI,
                                      PixelFormats.Gray32Float,
                                      bp,
                                      plane,
                                      xsize * stridefactor);
        }
    }
}
